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ABSTRACT 



A bulk turbulence-closure mixed layer model is generalized to allow prediction of very 
deep polar sea mixing. The model includes unsteady three-component mrbulent kinetic 
energy budgets. In addition to terms for shear production, pressure redistribution, and 
dissipation, special attention is devoted to realistic treatment of thermobaric enhancement 
of buoyancy flux and to Coriolis effects on turbulence. The model is initialized and 
verified with CTD data taken by R/V Valdivia in the Greenland Sea during winter 1993- 
1994. Model simulations show (i) mixed layer deepening is significantly enhanced when 
the thermal expansion coefficient’s increase with pressure is included; (ii) entrainment rate 
is sensitive to the direction of wind stress because of Coriolis; and (iii) the predicted 
mixed layer depth evolution agrees qualitatively with the observations. Results 
demonstrate the importance of water column initial conditions, accurate representation of 
strong surface cooling events, and inclusion of the thermobaric effect on buoyancy, to 
determine the depth of mixing and ultimately the heat and salt flux into the deep ocean. 
Since coupling of the ocean to the atmosphere through deep mixed layers in polar regions 
is fundamental to our climate system, it is important that regional and global models be 
developed that incorporate realistic representation of this coupling. 
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I. INTRODUCTION 



A. BACKGROUND 

All interaction of the atmosphere with the ocean must occur where the two fluids are 
in contact - that is, in the upper, fully turbulent layer of the ocean known as the mixed layer. 
This mixed layer is typically more buoyant than the water below, so that changes in the 
atmosphere are to some extent isolated from the deep ocean by the more buoyant layer. 
Exceptions to this general structure are in the polar ocean regions, where mixing may occur 
to great depths, or even to the bottom. This coupling of the ocean to the atmosphere through 
deep mixed layers in polar regions is fundamental to the initiation of the global thermohaline 
conveyor belt and to the earth’s climate system. 

Mixed layer modeling has historically used a linear dependence of density upon 
temperature, so that the thermal expansion coefficient is treated as a constant. This 
approximation is valid over much of the ocean, where mixed layers eire typically shallow, but 
breaks down in the ocean regions where there are large excursions in mixed-layer depth 
and/or mixed-layer temperature. The equation of state for seawater shows that the 
dependence of density and buoyancy upon temperature is not simply linear; the thermal 
expansion coefficient varies with pressure, and does so most at temperatures near freezing 
(Garwood, 1991). The thermobaric effect on the stability of a water column is most 
significant then where nearby parcels - very cold and at depth - have strong contrast in their 
potential temperatures. This is a characteristic state common to deep mixed layers in the 
polar seas, where traditional mixed layer modeling has failed to predict realistic mixing and 
deep water formation. 
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Another common approximation that may not be valid for deep mixed layers is the 
neglect of planetary rotation’s effects on turbulence. In more shallow mixed layers, this 
neglect is reasonable because the scale of the turbulence is much smaller than the distance 
over which the effects of rotation can be felt, the Rossby radius. However, in polar regions 
the Rossby radius is much reduced (~1 km). The scale of turbulent eddies is equivalent to 
the depth of mixing, so if deep mixing occurs (~1 km), the eddies are large enough to be 
affected by planetary rotation, and the effects need to be included in models of deep 
entrainment. Furthermore, the vector planetary rotation (o = ) needs to be 

considered because of the three-dimensionality of turbulence. 

B. OBJECTIVE 

Deep mixing is controlled by factors more complex than those of traditional mixed 
layer physics; the proper representation of non-linear buoyancy effects, energy-converting 
rotational effects, non-steady atmospheric forcing and water column initial conditions is 
critical to calculating the fluxes of heat, salinity, and tracers into the deep ocean. In this 
research a model was developed and tested that expands on traditional mixed layer modeling 
by adding the effects of thermobaricity and planetary rotation, that allows initialization with 
hydrographic data, and that forces the mixing with realistic meteorological time series of 
surface heat flux and wind stress. The goal is not just to compare the final mixed layer 
depth to that observed, but to see how the observed event is simulated in terms of relative 
contributions of energy and evolution in time. 

Chapter II summarizes the state of the art leading up to this research. The model 
equations are developed in Chapter III, and the numerical approach is discussed. Chapter IV 



2 



presents the oceanographic data used to initialize and the meteorological data used to force 
the model. Chapter V shows model results for hypothetical cases with constant forcing, 
followed by results with time dependent forcing. Chapter VI summarizes the results of the 
research and makes conclusions and recommendations for further study. 
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II. MIXED LAYER THEORY 



A. GENERAL MIXED LAYER AND DEEP WATER FORMATION THEORY 

Mixed layer theory concerning normal mixing depths (<100 m) has been well 
developed and is summarized by Zilitinkevich et al. (1979), Garwood (1979), and more 
recently by Large et al. (1994). However, none of these has direct application to very deep 
mixing (order 1 km and greater). 

Early speculations that deep water is formed not just by drainage off continental 
shelves, but by convection in the open ocean, were made by Nansen ( 1 9 1 2) for the Arctic and 
by Wiist (1928) for the Antarctic; it was not until much later that there was enough evidence 
to prove that open ocean deep convection does occur. 

Brown and Beardsley (1978) suggested that overturning during winter storms could 
be the primary source of mixing that forms Gulf of Maine Intermediate Water. They were 
cautious at the time, however, because they felt their data did not completely prove this 
result. 

Gordon (1978) reported the discovery of a very localized region (~13 km) of 
convection in the Weddell Sea. Killworth (1979) showed that the entire region of the 
Weddell Sea gyre was susceptible to overturning, and that reduced stratification within 
cyclonic eddies could be the mechanism that “selects” such narrow areas of convection 
within the larger gyre. 
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B. THERMOBARICITY 



Gill (1973), considering the mixing of water flowing off the continental shelf into the 
Weddell Sea, noted the importance of the variation of the thermal expansion coefficient with 
depth in cold water. He included this effect in his calculations and determined that there 
could be parcels of water with potential temperature - salinity ( 9 - s) properties such that, 
if displaced a finite distance downward, the parcels would be unstable, their relative density 
to the surrounding water continuing to increase as they descended, thus accelerating the 
parcels downward. 

McDougall (1987) coined the term thermobaricity for the increase in seawater’s 

thermal expansion coefficient with pressure (p ), or the ^ term in the equation of 

op oQ 

state, investigated its effect on mixing across neutral surfaces, and discovered that although 
both upward and downward vertical velocities across the neutral surface are possible, they 
are usually downward. 

Garwood (1991) studied thermobaricity’ s effect on mixing across a non-neutral 
surface, the interface at the bottom of the mixed layer. Although previously neglected, the 
buoyancy-flux enhancement by thermobaricity was found to be an important source of 
vertical turbulent kinetic energy, leading to increased entrainment. 

Garwood et al. (1994) noted that in addition to enhancing entrainment rates, 
thermobaricity might assist in bringing surface water to depth by causing water column 
instabilities. A parcel-type instability would occur when a parcel from the cold and fresh 
turbulent layer overshoots the interface, and continues to accelerate downward because 
pressure compresses it more than the surrounding water is compressed. A layer-type 



6 



instability would occur as an interface advected downward becomes unstable due to the 
thermobaric effect. 



C. ROTATION IN ATMOSPHERIC AND OCEANIC BOUNDARY LAYERS 

Rossby and Montgomery (1935) were the first to suggest that geophysical boundary 

layers should scale vertically with , where u, is the friction velocity and f=2Q^ is the 

vertical component of planetary rotation. Hence a boundary layer “Rossby number” is 
u 

— - 1 . 
hf 

Gascard (1973) showed that in a homogeneous water column, large amplitude 
internal waves generated by wind surges in shallow regions can propagate upward because 
of planetary rotation, bringing energy with them, and suggested they may couple with q and 
be involved in mixing and deep water formation. 

Wyngaard et al. (1974 ) were the first to evaluate turbulence terms containing ; 
their study was of the convective atmospheric boundary layer, and found ’s effect to be 
small. However, turning to ocean turbulence, Garwood (1977) found a Rossby number 
dependence for oceanic mixed layers in a bulk turbulence closure model. Garwood et al. 
(1985) went on to show that planetary rotation acts to exchange turbulent kinetic energy 
between horizontal and vertical components in a steady-state mixed layer. Vertical turbulent 
energy and thus mixed layer depth are increased under easterly winds. Garwood (1991) 
extended this hypothesis to the more general case of deepening mixed layers. 
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III. MODEL 



A. DERIVING THE EQUATIONS FOR TURBULENT KINETIC ENERGY 



Representing the Navier-Stokes equations of motion as 



dt 



^ dx^ 



1 dp 
P dx, 



- 2e. ..Q u. 

ijk j k 






d^u. 
1 

dx. 



( 1 ) 



the quantity turbulent kinetic energy is formed by decomposing the variables into their mean 
and turbulent parts, for example, representing one of the three velocity components as 

( 2 ) 

where the turbulent part is not assumed to be small compared to the mean. Here i = 1,2,3 
corresponds to the Cartesian coordinate directions of east, north, and up; 
[QJ =Q[0 COSO) s in 0 ) ]; molecular viscosity is v; gravitational acceleration is g; and 
<t> is latitude. The u-component kinetic energy per unit mass contained in the turbulence will 
then be 



Since 



_a_ 

dt 




du du 



( 4 ) 



and the right-hand-side contains only quantities predicted by the equations of motion, kinetic 
energy can be calculated. Again generalizing to the three components, the last term can be 
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written as 



— du du _ du^ _ 1 ap 

^ dt ^ ^ dx. ^ ^ dx. ^ po dx. 



u . Q. 

dx, Po 



( 5 ) 



after Reynolds’ decomposing, averaging, and making the Boussinesq approximation on 
equation (1). The first term on the right-hand-side of equation (4) is formed by performing 
the multiplication and decomposition first, and then averaging, resulting in 



du. 



du . 
1 



du; 



u^du^ du. 



■ Uv u. - u,- u. 



dui 



dt dx^ dx. dx. dx. dx. 



1 — dp 
Po ' dx, 



1 , dp' , — 

- Uj^ -r^— + V ( U 

Po dx, 



d^ u . d^ u ■ 

+ u; ^ ) 






dx‘ 



~^^xjk S ( u. p - u: p' ) (6) 



Subtracting, we have 



I du^^ du; du, 

”2 at 



au; 



= - u. Uv 



1 , dp 



^ dx. dx. dx^ Po dx. 



a'u; 



+ V u . u 



. s-,- 2 e . .. Q . Uy Uju “ 5 — Uy p' 

dx. Po 



( 7 ) 



The first and third terms on the right-hand-side are simplified by using the reverse product 



rule for differentiation and the incompressible form of the continuity equation. Then 



and 



— , du; I s 



u, u, 






^ ^ dx . 2 dx . 

j j 

The TKE equation for the system is 



( 9 ) 



at 



d u . 



dx. 



- 2 a:u;^ 



2 dUi 2d, — , 

— P'~ ^ ( UiP' ) +^^Ui 

p dx p dx 



d-u' 



dxj 



4 e 



ijk 



Q u; u; -2 6 _j3^u;p' 



(10) 



Changing notation so that [ u J = [ u, v, w ] , [xj = [ x, y, z ] and assuming 



^ 

horizontal homogeneity so that — ( ) 




for tvirbulent kinetic energy: 



gives the component equations 



du-^ 

dt 



2 w w‘ 



du 

dz 



az 



2 du 
+ — P' 

Po 



- ^Q^U' w' + 4Q^u'v' 



+ V ( li 



d^ u' 



dx^ 



d^u' 

dy^ 



+ w 



dz^ 



( 11 ) 



at 



- 2 



dv 



V' w — — 
d z 



- A ( „- v'M 

d z 



2 dv^ 
+ — P' 

Po 



AQ^u' V' 



+ V 



d^ V' d^ V' 

V' + V' 

dx^ dy^ 



d"v' 

dz^ 



( 12 ) 
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dw^ 

dt 



- 2 —p' w 

Po 




+ 4 n.W w' + V { w‘ 



w 

dx^ 



d^w 

dy^ 



d^w’ 

dz^ 



( 13 ) 



where the terms on the right-hand-side are shear production, turbulent transport, pressure 
redistribution, pressure transport, buoyant production/damping, Coriolis redistribution, x-z 
rotational redistribution, and viscous dissipation; and the term on the left-hand-side is called 
the storage term. 

B. INTEGRATING FOR BULK TKE 

Integrating over the depth of the mixed layer plus entrainment zone will give the 
bulk TKE for the layer. The storage term becomes 



where < > denotes the vertically integrated value. Turbulent transport and pressure 

transport terms integrate simply to zero, assuming no turbulent energy or pressure 
perturbations are transported across the boundaries. (Entrainment does not cross the 
interface at z = -h; rather, it changes the depth of h.) 



Integration of the buoyant production/damping term and the x-z rotational redistribution term 
requires consideration of the approximately-linear flux profiles of temperature, salinity, and 
momenttim within the layer, and the “jump condition” at the bottom of the layer. For any 




(14) 




(15) 




(16) 
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conserved quantity C, the linear flux profile has a shape such as Figure (1) and the jump 
condition is 



C' W' \ .= ~ W A C 
-n e 



(17) 



This introduces a new variable w^, the entrainment velocity. In the absence of upwelling or 
downwelling, 



dt 



(18) 



Alternatively, if non-zero mean vertical velocity is included, then 

dh 



w - w \ ^ 
dt ^ 



The temperature flux at the surface is 



T’ W I 






(19) 



( 20 ) 



where is the rate of cooling of the surface, and is the specific heat of seawater at 
constant pressure. The surface salinity flux is 



S' w' I 0 = S Vol 



( 21 ) 



where Vol is the rate of addition or removal of water from the surface, the net result of 
precipitation, evaporation, melting, and freezing, expressed in units of length per time. The 
equations for the fluxes as functions of z then become 



and 



T W (z) 



Qq 

PoS 




-w AT ( — ) 
® h 



( 22 ) 



S'W'(z) =SVol{l*—) -wAS{— ) (23) 

h ^ h 

Rather than treating the thermal expansion coefficient as a constant, it is allowed to vary 
linearly with depth. 

a (z) = a^- a^z (24) 
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Figure 1 . Flux profiles shown are for momentum, but would have similar shapes (piecewise 
linear) for any conserved quantity. The upper plot is for wind direction west to east, and the 
lower plot is for wind east to west. 
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Then the integration of the buoyant production and x-z rotational redistribution terms reduces 
to evaluating definite integrals of polynomials, so 

f° -2 — p'w dz = 2 f° b' w dz 

J-h Pq J-ft 



= 2 f° a (z) T' w' dz - 2 f° H S' w' dz 
J-h J-tt 



0 , _ 



= gh[ ^ ( a„ + 



) + 3 s voi 



a, h 



- ghw [ AT( a + ) - 3 AS ] 



1.5 



(25) 



and 



f° 4Q. u' w'dz = 2Q,h ( — + w At; ) (26) 

Po 

By scaling, the second term on the right-hand-side of equation (26) is small and is neglected 
hereafter. 

The integral of the shear production term has contributions only from the regions 
where there is velocity shear, near the surface and just below the bottom of the fully turbulent 
layer, as shown in Figure (2). Measurements have shown the velocity shear just below the 
surface to be proportional to the friction velocity, so that 

5u = m^u^sinip (27) 



where <p is the angle from north clockwise to the v^nd vector, , and 



jn. 



is the constant of proportionality. The momentum fluxes at the surface and within the 
entrainment zone are approximated by the piecewise linear functions in Figure (1). 
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Figure 2. To model the production of turbulence by shear in the mean flow, values of 
vertical shear at the ocean’s surface and in the entrainment zone are estimated rather than 
assumed infinite. 
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The surface momentum flux w w | and the shear near the surface | are always 
opposite in sign, so that their product is always negative. Using the relationship 

^ = uj-sin<p (28) 

Po 

the integral becomes 



- 2 /° 
J -h 



— — T du 
u ' w' dz 

-6 0 Z 



- 2 I u' w' dz 
J-6z dz 



o r-h ~ — r du , 
-2 1 u' w' dz 
J-h-b dz 



= -2 u' w' L 5l2 -2 u'w' I 6 Au 

U -ft- — 



2 in. X ^ _ 

( — ) ^ + w(Au) 



(29) 



Considering the Coriolis redistribution, the Reynolds’ stress is parameterized as 
proportional to the horizontal shear of the mean flow, which under horizontal homogeneity 
is zero. This means that because both forcing and response are the same everywhere, there 
can be no net transport of u turbulence in the v direction, and no net transport of v turbulence 
in the u direction. 

4Q^ f° w V dz “ ( ) = 0 (30) 

'‘J-h oy ox 

Using the “return to isotropy” assumption of Rotta (1951), the pressure redistribution 
term becomes 

= 2m. ( <E> - 3<u^> ) <E>'^ (31) 

Pgj-fl dX 



where <E> = ( <w^> + <v'^> + <w^> ) and is a constant of proportionality. 
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Applying the product and chain rules for differentiation, the dissipation integral has 



terms of two forms. 





du' du' 



du' du' 




V 



+ 



+ 



dx dx 



dy dy 



+ V 




(32) 



The last term, viscous diffusion of turbulent energy, is negligible compared to both buoyant 
and shear production, or to turbulent transport. The remaining terms on the right-hand-side 
of equation (32) are u-component viscous diffusion and represent all of the dissipation of u 
TKE. Total dissipation is assumed to be a function of the net turbulent energy only, and is 
thus; 



where is a constant of proportionality and < > denotes the mixed-layer average value. 

Since viscous dissipation occurs at the smallest scales, where kinetic energy is converted to 
heat, total dissipation is also assumed to be isotropic. The corresponding terms of each 
component equation would then be one third of the total dissipation. 

Using the integral results of equations (14), (25), (26), (29), (31), and (33), the bulk 
equations for the components of turbulent kinetic energy are: 



3 



D = 2m, <E> ^ 



(33) 
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2 

+ 2m2 { <E> - 3<u'^> )<E>^ 



(34) 



- -m,<E> 



2 

2 



dh<v'^> 

at 



2m^ T ^ — 

i (-^)^ + wi/\v) 

Po 



2 

+ 2/n2 ( “ 3<v'^> )<E>^ 






2 

2 



(35) 



dh< w'^> 
at 



g-h [ 





+ P S Vol 



ah 

- ghw^[ AT( 3g+ y-g. ) - 3 AS ] 



- 2Q.h — 
Po 



2 

+ 2m^ ( <E> - 3<w'^> )<E>^ 






(36) 



C. TEMPERATURE AND SALINITY EQUATIONS 

Conservation principles, plus horizontal homogeneity, alloAV the potential temperature 
and salinity equations to be written as 

ae 
at 

and 

at 

Integrating over the mixed layer plus entrainment zone, and applying the jump condition, 
equation (17), the equations become 



- ae 

w—. 

dz 



a©' w' 

a z 



-ds 

w — 

a z 



as' 



a z 



(37) 



(38) 
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and 



ae ^ ^ gp 

dt h 



- w^AQ ) 



(39) 



4^ = 4 ( 5Voi - W AS ) 

d t n 



(40) 



D. BULK MOMENTUM EQUATIONS 



Equation (1), after decomposing, averaging, making the Boussinesq approximation, 
assuming horizontal homogeneity, and neglecting viscous effects on the mean flow, gives 



the component equations for mean flow: 



d u 
It 




— 4 ^ - 2 Q^w + 2 n.v 

p^ dx 



du' w' 
d z 



(41) 



dv — dv 1 dp — dv* w' 

— = - ^ - 20 U - r 

dt dz Pq dy ^ dz 



(42) 



dw 

Tt 



— 1 dp 

dz Pq 






P-Pn 



+ 20^ 



dw‘ w' 
dz 



The hydrostatic approximation is made, so that equation (43) becomes 



3p 

dz 



= -pg 



and upwelling is neglected here. 



(43) 



(44) 



Although the terms were retained in the turbulence equations because 
w’ - o ( u ' ) , they can be dropped here because w « 0{u,v) . Also if the geostrophic and 



wind-driven flows are separated, equations (41) and (42) become 



dt 




+ 2Q3V 



dw w 
dz 



(45) 



and 
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dV 

dt 



( 46 ) 



where u and v represent the wind-driven flow averaged over some short period of time. 
Integrating equations (45) and (46) yields 



and 



d{h<U>) 

dt 



= 2Q^h<V> 




d (h<V>) 
d t 



= - 2Q^h<U> 



(47) 



(48) 



E. ENTRAINMENT VELOCITY 

Entrainment velocity must now be parameterized to close the system of equations; 
the local TKE budget at the bottom of the mixed layer, where the entrainment is occurring, 
can give a scale estimate. The sum of equations (1 1), (12), and (13) retains contributions to 
the storage term from shear production, turbulent transport, buoyant damping, and 
dissipation: 



d i u'^ + v'^ + w’^ ) 
d t 
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Po 



(49) 
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where e = 



V ( U ' 
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dx^ 



+ u' 



dy^ 



+ L7' 



a u' 
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a^ ^ 



dx^ 



V' V' 

V' + V' 

dy^ dz^ 



w* 



a^ w' 



w' w' . 

w' ^ + w' ) 



(50) 



dx^ dy^ dz^ 

Near z = -h both shear production and dissipation are small compared to the turbulence 

transported from above and to the amount being removed by buoyant flux, so that the balance 
is largely among storage, turbulent transport and buoyant damping. 



d { u'^ + v'^ 



w'^ ) 



dt 



^ [ W' {E + 2-^) ] - 2 ^P' w 
p„ p„ 



(51) 



Each of these terms can be scaled on bulk variables, and the time scale for changes in 
turbulent energy at the bottom of the mixed layer can be thought of as the time required for 
a change in energy content to be turbulently diffused from the surface, which is posed here 
to depend upon h and <w^> . 

(52) 



+ V'^ + W'^ ) 



d t 



:£> V 



-^ [ w [E * 2^) ] = 

dz Pq h 



<w'^> 



(53) 



2 ^p-w’ 



■ w gh ( aAT - & AS ) 



(54) 



Combining and solving for entrainment velocity give 



W = 
e 



<E> y <w’^> 
<E> + gh ( aAT - &AS) 



(55) 
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F. 



NUMERICAL APPROACH 



The system of eight equations, summarized as Table (1), was cast as a MATLAB 
function, and the function was then called using the MATLAB utility ODE45, a Runge-Kutta 
ordinary differential equation solver that has automatic step-size control. The tolerance for 
the solver was relaxed from its default value of 10'^ to 10"^ for faster convergence. 
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Table 1 . Equations Used in the Model 
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IV. OBSERVATIONAL DATA 



A. IMPORTANCE OF THE GREENLAND SEA REGION 

As summarized by Visbeck et al. ( 1 995), the Greenland Sea has long been considered 
to have characteristics favorable for deep water formation. Nansen, in 1906, observed only 
small differences in properties of Greenland Sea surface and deep waters, and proposed that 
deep water formation occurred there. Similar (but more modem) observations were made 
by Rudels et al. (1989), when mixing seemed to have recently extended to 1500 m. A 
s}mthesis of the physical explanations provided by Visbeck et al. (1995) and Schott et al. 
(1993) for how the region becomes favorable to deep convection includes the influences of 
the previous seasons as well as local forcing. Bourke et al. (1992) showed that during the 
ice-free months of the year, upper water conditions within the basin can be variable due 
largely to advection by the Jan Mayen Current. In early winter, as the marginal ice zone 
propagates rapidly eastward, brine release and entrainment act to increase the salinity of the 
surface waters, beginning to erode the buoyancy of the mixed layer. Deepening during this 
period may extend far enough into the warm Atlantic Intermediate Water, so that when 
subsequent cooling and stirring occur, warming by entrainment will dominate over surface 
cooling; the formation of ice is prevented and an ice-free area, surrounded on three sides by 
ice but open to the north, known as the Nordbukta, forms. The Nordbukta, situated over the 
Greenland Sea gyre, is then preconditioned for deep convection. Hence, the data discussed 
below, provided by Gascard and Lherminier (personal communication), were very well 
suited for use as a test case for the model. 
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B. 



OCEANOGRAPHIC DATA 



Two station soundings taken by the RA^ Valdivia at 74.5 °N, 2.5° W, during 
Februaiy-March 1994 were examined. At this station, the ocean depth is more than 3600 m. 
For each sounding, the pressure, depth, in situ temperature, potential temperature, salinity, 
potential density referenced to the surface, potential density referenced to 2000 m, sound 
velocity, and dynamic height were available. In addition to the oceanographic data, each 
sounding had recorded the station number, date, hour, precise latitude and longitude, depth 
of sounding, sea surface temperature, sea surface salinity, wind speed and direction, and air 
temperature. Figure (3) shows February 16th’s vertical profiles of potential temperature, 
salinity, and potential density referenced to the surface. A 200-m deep mixed layer of cold 
fresh water overlies more saline and warmer water; the buoyancy effects of temperature and 
salinity are in opposition, but salinity dominates. The 16 February profiles show that the 
magnitudes of the temperature and salinity gradients at the bottom of the mixed layer are 
very large, exceeding 0.1 C/m and 0.01 psu/m, respectively. Figure (4) shows profiles of the 
same quantities for March 19th, 31 days later; no intermediate soundings were available. 
The mixed layer is much deeper (600 m), and the entrainment of warmer more saline water 
has diminished the contrast between the two water masses; the buoyancy jump has been 
eroded by the mixed layer deepening. 

Below 600 m the two soundings look very similar, especially if the March sounding 
is offset upward by about 30-40 m. This suggests a slow but steady downwelling speed of 
about 1 m/day, which is small compared to the average entrainment speed of about 1 3 m/day. 
Comparison of heat and salinity content between the tw'o soundings gives an average loss of 
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Figure 3. Sounding of 16 February 1994 shows a 200 m deep mixed layer overlying warmer, 
more saline water. 
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Figure 4. Sounding of 19 March 1994 shows mixed layer has deepened to 600 m. 
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130 W/m" and no surface flux in salinity. The similarity of profiles at depth as well as the 
apparent conservation of salinity gives some indication that horizontal advection is small, 
and that one-dimensional modeling should be useful for this case. As will be seen in Chapter 
V, the calculated Eckman transport for the layer suggests one-dimensional modeling as well. 

C. METEOROLOGICAL DATA 

Meteorological data from the European Center for Medium Range Weather Forecasts 
(ECMWF) were available for February 2nd through May 3 1st. Parameters were six-hourly 
values of surface temperature, 2 m temperature, 10 m u-component of wind, 10 m v- 
component of wind, surface sensible heat flux, and surface latent heat flux. The values were 
interpolated from nearest-neighbor 2° grid points in the ECMWF analyses. Figure (5) shows 
the wind speed and associated wind stress for the period February 16th through March 19th, 
and Figure (6) shows the surface latent and sensible heat fluxes for the same period. The 
mean heat loss for the period indicated by these meteorological estimates is 1 82 W/m^, 
significantly larger than the actual average heat loss measured in the water column. This 
difference may be explained by lack of inclusion of net radiation or by error attributable to 
coarse resolution in the ECMWF product. Methods of correcting for this difference are 
discussed in Chapter V. Three events stand out as periods of large surface forcing, or storms; 
the latest of the three is the largest. It would be interesting to determine whether the mixed 
layer depth was relatively stable until this final large storm mixed the layer down to 600 m, 
or whether the deepening occurred gradually. 
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Julian Day 



Figure 5. Meteorological data from ECMWF shows that wind speed and wind stress are 
highly variable over the period. 



Watts/m^2 



SSHF SLHF 





SSHF + SLHF 




Figure 6. i hree episodes of large surface heat loss stand out in the sum of surface sensible 
and latent heat fluxes, from ECMWF. 
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V. MODEL RESULTS 



A. EFFECT OF THERMOBARICITY 

The general behavior of the model was tested, comparing results in cases where 
thermobaricity is included to cases where it is neglected. These cases demonstrate that the 
increase in thermal expansion with pressure (depth) does increase the depth of mixing. 
Figures (7) and (8) show the deepening and the TKE components for a 25-day period with 
constant forcing of southerly winds (0.8217 dynes/cm^ wind stress) and moderate cooling 
(100 W/m^). Including the thermobaric effect increased deepening by 50%. Figures (9) and 
(10) show the deepening and the TKE components for a 25-day period with the same forcing 
except stronger cooling (400 W/m^); including the thermobaric effect increased deepening 
by 33%. The point where the vertical turbulent energy exceeds the wind-driven turbulent 
energy (v-component in this case) can be interpreted as a shift from forced-convection 
dominant to free-convection dominant. The shift of regimes occurs earlier when 
thermobaricity is included. 

Parcels above and below a typical Arctic region mixed-layer interface have a smaller 
buoyancy contrast with o than they have with thermobaricity neglected; they also have 
a smaller buoyancy contrast at depth than they would at shallower mixed layer depths. This 
diminished buoyancy jump will act to increase w^, increasing the deepening rate. Further, 
the buoyancy flux at the bottom of the layer is proportional to both the entrainment velocity 
and the buoyancy jump; entrainment buoyancy flux is increased by both of these factors. 
Another process that deepens more due to thermobaricity is surface cooling; heat taken 





seconds x10® 

Figure 7. With 100 W/m‘ of cooling, the model deepens the mixed layer 50% more when 
thermobaricity is included in the calculation (a, * O). Wind stresses Taux and Tauy are in 
dynes/cm"; heating QO is in W/m‘; and a, is in (cm 
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Figure 8. Velocity variances ( u , for example) represent the turbulent energy present in 
the eddy field of each component. Comparison of the components in the plots above shows 
that the increase in deepening for the case with thermobaricity is due to increased turbulent 
energy in the vertical component. 
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Figure 9. With 400 W/m* of cooling, the model deepens the mixed layer 33% more when 
thermobaricity is included in the calculation (a. O). Wind stresses Taux and Tauy are in 

dynes/cm‘; heating QO is in W/m-; and a, is in (cm °C)''. 
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Figure 10. Turbulent energy components, represented by velocity variances, show that the 
increase in deepening is due to more vertical turbulent energy. The shift from forced 
convection to free convection occurs earlier with thermobaricity included. 



37 



out at the surface will have a greater density effect at depth than at shallower depths. 
Thermobaricity increased the deepening rate even for cases with as little as 1 0 W/m‘ of 
cooling. 

B. EFFECT OF COOLING 

To examine the sensitivity to seasonal cooling, the model was integrated for a period 
of 100 days from the initial profile of 16 February, using forcing corresponding to the means 
of u-wind and v-wind and two different rates of cooling. In the first case, the observed 
average heat loss of 130 W/m* was used; in the second case, twice that. Figure (1 1) shows 
that doubling the rate of cooling has the effect of doubling the deepening rate; all other initial 
conditions and forcing are the same, and the heat is simply removed more quickly. Figure 
(12) shows that free convection dominates much sooner with the stronger cooling. It is also 
interesting to see the response of the horizontal turbulent energy components to the dramatic 
increase of vertical turbulent energy. They increase as a result of pressure redistribution; the 
turbulence is attempting to become more isotropic. 

An interesting question is whether it was possible for the water column as measured 
on 16 February to be mixed to the bottom before summer ice melt and warming shut off deep 
mixing. Looking at Figure (11) again with this question in mind, it seems that there could 
not be enough time left even with very strong cooling for mixing to the bottom to have taken 
place in 1994. 

C. EFFECT OF ROTATION 

The direction of the wind has a potential effect on deepening because it can oppose 
or coincide with the direction of planetary rotation. To examine this effect, two 100-day 
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Figure 1 1. Doubling the cooling rate doubles the rate of deepening. 
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Figure 12. ihe vertical component of turbulent energy is dominant much sooner when 
cooling is 260 W/m" vice 130 W/nr. 
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cases were simulated, one with large westerly wind stress (+2 dynes/cm^) and one with large 
easterly wind stress (-2 dynes/cm"). Figure (13) shows that the westerly case deepened to 
1289 m, while the easterly case deepened further, to 1468 m. In Figure (14), while the total 
turbulent energy for the two cases is about the same, the amount of energy that is in the 
vertical TKE is strikingly different. With easterly wind, u turbulent energy is converted to 
w turbulent energy (Figure (15) ); with westerly wind, w turbulent energy is diminished by 
conversion to u turbulent energy (Figure (16) ). The lower total energy for the deeper case 
is explained by the fact that there has been more buoyant damping (sometimes called buoyant 
consumption) associated with the increased entrainment. 

D. CONSTANT AVERAGE FORCING RESULTS COMPARED TO DATA 

When the model is integrated for 31 days with the mean wind components and 
observed average cooling, a mixed layer depth is calculated that is very close to that observed 
(627 m). However, averaging the vector wind components under-represents the amount of 
stirring; in order to make a more valid comparison, mean wind speed was computed and a 
typical wind direction chosen. Components of wind stress were then calculated from these, 
and the depth of mixing after integrating for 3 1 days was 733 m (Figure (1 7) ). 

E. TIME DEPENDENT FORCING 

The model was solved for the 3 1 days between the times of the CTD casts of 1 6 
February and 19 March by linearly interpolating the winds and heat fluxes for times between 
the six hourly estimates from ECMWF. A final depth of 802 m was reached, and Figure (1 8) 
shows that the periods of more rapid deepening correspond to the more intense forcing 
events. The deepening of the mixed layer did not occur in one large surge during the third 
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Figure 13. Westerly winds (from the west, upper plot) are not able to deepen the mixed layer 
as much as easterly winds (lower plot). Wind stress is in dynes/cm*; heating is in W/m". 
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Figure 14. Total turbulent energy is nearly the same for the two cases, but the amount that 
is in the venical eddies is different. Wind from the east (lower plot) puts more energy into 
the vertical component. 
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Figure 15. With forcing by easterly winds, the contribution from rotation is negative for u 
TKE and positive for w TKE: w receives energy from u. 
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Figure 16. With forcing by westerly winds, the contribution from rotation is positive for u 
TKE and negative for w TKE: w gives up energy to u. 
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Figure 17. The upper plot shows deepening over 31 days with forcing computed from mean 
wind components. The lower plot shows deepening over 3 1 days with forcing computed from 
mean wind speed and a typical direction. Wind stress is in dynes/cm'; heating is in W/m'. 
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Figure 18. Mixed layer deepening is calculated with time-dependent forcing from 
meteorological data, and is plotted above graphs of wind speed and heat flux. The periods 
of more rapid deepening correspond to the periods of more intense forcing. In this model run 
the heat fluxes were corrected by scaling. 
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and strongest storm, but happened gradually. Even though the final forcing event was 
significantly larger than the other two, the layer was more resistant to deepening because it 
was already deeper. For a deeper layer, more energy is required for the same amount of 
entrainment simply because there is more water that must be kept turbulent. 

Although the magnitude of heat removal from the ocean as indicated by the 
meteorological data was too high, there is a high degree of confidence in the qualitative 
correctness of that data; the storms’ occurances, their relative magnitudes, and the phase 
relationships of wind peaks to heat flux maxima are taken to be well represented. In other 
words, the shape of the surface heat flux time series should be retained even though its 
magnitude should be adjusted. Offsetting the heat flux by uniformly subtracting 52 W/m^ 
results in the correct total heat loss from the ocean, but results in periods of net warming of 
the ocean by the atmosphere in winter which are thought not to be realistic. Uniform scaling 
( multiplying by ) did not present this disadvantage, and so was preferred. Another 

option of first adding 50 W/m’ as an approximate additional heat loss from the ocean due to 
net long wave radiation, then multiplying by , was investigated later, but the 

results were so similar as to be almost indistinguishable. 

Since momentum budgets were calculated in the model, they provided another check 
on the reasonableness of one-dimensional treatment of this problem. Integrating the mean 
wind-driven U and V of the mixed layer over each time interval, the mean wind-driven 
transport is computed. The progressive vector diagram for this transport is shown as Figure 
(19), and indicates the water column would move ~5 km to the west, 90 degrees to the right 
of the mean wind direction, which is almost due south. 
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Figure 19. The time integral of the wind-driven velocities U and V of the mixed layer is 
plotted as a progressive vector diagram, showing the inertial effects and the mean transport 
of just over 5 km almost due west. 
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The components of turbulent kinetic energy are represented by , and 

and are plotted in Figure (20). What is striking in that figure is that v TKE was very 
dominant while w TKE accounted for only 10 to 20% of the total TKE during storm events; 
the turbulence was very anisotropic. 

Figure (2 1 ) shows the contributing terms to u TKE. Shear production at the interface 
is always small, and pressure redistribution is the same order as shear production at the 
surface. Figures (22) and (23) are the corresponding plots for terms contributing to the v and 
w TKE budgets. In Figure (22), shear production at the surface is large, and v is losing even 
more energy through pressure redistribution than through dissipation. Figure (23) shows that 
the main source of vertical turbulent energy was pressure redistribution of wind energy'. Also 
interesting is that the entrainment buoyancy flux is often as large as the surface buoyancy 
flux, and is plainly not a fixed fraction thereof* The integrated effect of rotation over the 
period is neutral, but rotation’s effect during any particular storm may not be. Figure (24) 
is an enlargement of one section of the time series; the influence of rotation during the first 
event was to enhance entrainment, and during the second to inhibit entrainment. These 
observations show that for the same magnitude of forcing, a different depth of mixing may 
result. The phase relationship between wind stirring (via pressure redistribution) and surface 
buoyancy flux varies for different events as well. Could this also be a factor in the 
effectiveness of a particular storm for deep mixing? 



* A common approximation for entrainment by free convection in the atmospheric boundary layer is that 
the entrainment bouyancy flux is opposite in sign and a fixed fraction (order .2) of the surface buoyancy flux 
(Stull, 1988). 
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Figure 20. The turbulent kinetic energy components are plotted as functions of time for the 
time-dependent forcing case with surface heat flux scaled. Both wind stirring and cooling 
are vigorous, but v TKE dominates over w TKE. 
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u TKE Budget 




rhe contributions to u TKE are plotted as functions of time. They are: shear 
the surface, shear production at the interface, pressure redistribution, rotational 



, dissipation, and storage. 
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V TKE Budget 




Figure 22. The contributions to v TKE are plotted as functions of time. They are: shear 
production at the surface, shear production at the interface, pressure redistribution, 
dissipation, and storage. 
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w TKE Budget 
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Figure 23. The contributions to w TKE are plotted as functions of time. They are; surface 
buoyancy flu.x, entrainment buoyancy flux, pressure redistribution, rotational redistribution, 
dissipation, and storage. 
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Figure 24. A section of Figure (23) is enlarged to show detail. 
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Figure (25) shows the various terms contributing to the total TKE of the mixed layer. 
The oscillatory behavior of the storage term is an interesting depiction of the mixed layer’s 
ability to alternately store and then dissipate energy in response to wind surges. 

Turning to the question of why the final mixed layer depth on 19 March is ~200 m 
too deep, two basic but related factors are felt to be responsible: the model does not retreat, 
and it may transport energy to the bottom of the layer too efficiently. After a period of low 
winds, when the wind stress suddenly increases, the ocean mixed layer requires time not only 
to re-stir the layer before entrainment recommences, but to transport energy surges 
downward before entrainment velocity can increase. Also, just after increased stirring, 
energy that is concentrated near the surface of the ocean gives more dissipation than the same 
amount of energy would if it were distributed more evenly throughout the layer. All of these 
effects can be corrected for by adjusting the model constants: the dissipation constant m„ the 
pressure redistribution constant m,, the shear production constant m-, and the entrainment 
efficiency constant m 4 . However, in this research, model constants were not adjusted to 
achieve quantitative agreement since such an adjustment would have meant tuning a very 
general model to match one specific observed case. 



56 



evs/evLuo 



Total TKE Budget 




Figure 25. The contributions to total TKE are plotted as functions of time. They are: surface 
buoyancy flux, entrainment buoyancy flux, surface shear production, dissipation, and storage. 
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VII. CONCLUSIONS 



A. SUMMARY 

Coupling of the ocean to the atmosphere through deep mixed layers in polar regions 
is fundamental to our climate system, and much of the energy exchange occurs episodically 
during storms. In order to understand the physics of the ocean mixed layer response to 
storms, a bulk turbulence-closure mixed layer model was developed and was used to 
simulate an observed case of deep polar sea mixing. 

The model included budgets for potential temperattire, salinity, horizontal 
momentum, and unsteady three-component turbulent kinetic energy. In addition to shear 
production, surface buoyancy flux, entrainment buoyancy flux, pressure redistribution, and 
dissipation, realistic treatment was also given to thermobaric enhancement of buoyancy 
flux and to Coriolis’ effects on turbulence. The model demonstrated the ability to 
qualitatively simulate events when initialized with actual CTD data and forced with a 
realistic meteorological time series, provided the heat fluxes were adjusted to give a correct 
total heat loss from the ocean. Further work toward tuning of the model constants rr\, mj, 
mj, and m 4 , and the inclusion of a shallowing mode have the prospect of making the model 
quantitatively useful. 

Several conclusions can be made about the mixing that occurred in the Greenland 
Sea during late winter 1994: 

• Horizontal wind-driven energy exceeded vertical free convective energy throughout 
the period; the turbulence was anisotropic, and wind stirring was dominant over 
cooling even during a period of strong cooling. 
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• In the vertical energy budget, pressure redistribution of horizontal energy (from the 
wind) was a larger source of energy than surface buoyancy flux; free convection 
was not dominant even when considering only. This underscores the 
importance of including wind unsteadiness. 

• The deepening response for a given forcing event depends not only upon the 
strength of the event, but also upon the initial mrbulent energy, the depth of the 
layer, the buoyancy jump at the bottom of the layer, and the stratification below 
the turbulent region. 

• The direction of the wind may be important not only during strong events, but 
between them as well, because it influences the “spinup” time for the next 
entraining event. 

• The phase of the surface buoyancy flux to shear production also differs among 
events; it is possible that this difference may influence the storm’s efficiency at 
deepening the mixed layer. 

This study demonstrates the importance of water column initial conditions, the 
accurate representation of strong surface cooling events, and inclusion of the thermobaric 
effect on buoyancy to determine the depth of mixing and ultimately the heat and salt flux 
into the deep ocean. 

B. RECOMMENDATIONS 

One of the primary motivations for the study of deep mixed layers is to be 
eventually able to parameterize mixed layer dynamics, in particular the dynamics of deep 
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mixed layers and subgrid scale penetrative convection. The requirement is that the effects 
of these processes be included in general circulation models despite the impracticability of 
modeling the actual physical processes within the general circulation model itself. This 
research provides a tool for studying various forcing events and their calculated effects, 
to conceptualize and test what the important variables in such a parameterization might be. 
To increase this model’s usefulness, several areas need further study: 

• Include a shallowing mode capability. Turbulence may be dissipated even without 
warming, so that the fully turbulent layer may at times be much shallower than the 
layer that is well-mixed in potential temperature and salinity. The additional time 
required to re-stir the layer is currently neglected. 

• Include a finite thickness entrainment zone. Stability within or across the zone may 
be important. 

• Compare deepening and energy component behavior to results from large eddy 
simulations of the same events. Specifically, find better estimates of model 
constants, and ensure the correct amount of anisotropy. This would be an important 
step to improving the quantitative accuracy of the model. 
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